Defect physics and electronic properties of Cu 3 PSe 4 from first principles 
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The p-type semiconductor Cu3PSe4 has recently been established to have a direct bandgap of 1.4 
eV and an optical absorption spectrum similar to GaAs [Applied Physics Letters, 99, 181903 (2011)], 
suggesting a possible application as a solar photovoltaic absorber. Here we calculate the thermody- 
namic stability, defect energies and concentrations, and several material properties of Cu3PSe4 using 
a wholly GGA+C/ method (the generalized gradient approximation of density functional theory with 
a Hubbard U term included for the Cu-d orbitals). We find that two low energy acceptor- type de- 
fects, the copper vacancy Vcu and the phosphorus-on-selenium antisite Ps e , establish the p-type 
behavior and likely prevent any n-type doping near thermal equilibrium. The GGA+!7 defect cal- 
culation method is shown to yield more accurate results than the more standard method of applying 
post-calculation GGA+(7-based bandgap corrections to strictly GGA defect calculations. 

PACS numbers: 61.72.J-, 71.20.Nr, 71. 15. Mb, 88.40.fh 



The growing family of multinary copper chalcogenides 
has been of great interest for solar photovoltaic applica- 
tions. In addition to the commonly used solar absorber 
CuIni_ a; Ga a ;Se2 (CIGS), materials that have raised inter- 
est include Cu 2 ZnSnS 4 , Cu 7 TlS 4 \ CuClSe 2 2 , CuBiS^, 
CuSbS2 3 , and Cu3BiS3 4 . Recently the p-type semicon- 
ductor Cu3PSe4 has been established 5 - to have a direct 
bandgap of E g = 1.4 eV, with a calculated absorption 
a > 5 x 10 4 cm -1 for wavelengths less than 630 nm. 
This bandgap lies in the optimal range for photovoltaic 
power output and warrants further investigation of the 
material. 

In addition to optical absorption, essential considera- 
tions for photovoltaic applications include ease of syn- 
thesis, conductivity, amenability to doping, and trap- 
assisted charge recombination. These quantities are 
largely controlled by the thermodynamic stability of the 
material with respect to competing phases and point 
defects. Here we perform a point defect analysis of 
Cu3PSe4 combining the +U Hubbard term for total en- 
ergy calculations with the correction methods described 
recently by Lany and Zunger— Several potential sub- 
stitutional donor defects are also considered. Further- 
more we examine bulk properties including the partial 
density of states (DOS), the dielectric tensor, and the 
highly asymmetric effective mass tensor. We compare 
our results to recent experiments^ and to a more stan- 
dard GGA calculation. We also compare our methods 
with the alternative electrostatic image correction proce- 
dure described by Freysoldt et al£. 

Defect formation energies are most often calculated 
using density functional theory (DFT) within the local 
density approximation (LDA) or within the generalized 
gradient approximation (GGA). However, recent statis- 
tical studie d 10 ' 11 on the accuracy of heat of formation 
calculations indicate that using GGA with an additional 
Hubbard U term for the occupation of transition metal 
d orbitals, the so-called GGA+L 7 method, will be more 
accurate than using standard GGA or LDA. The U value 
is held constant for each type of transition metal atom 



throughout the analysis, including calculations of the en- 
ergies of the transition metal elements themselves. These 
studies also suggest that one should add a statistically de- 
termined correction value to the total energy of each pure 
element before calculating the heat of formation AH of 
a compound. To obtain the most accurate heat of for- 
mation energies for both compounds and defects, we use 
GGA+U— and apply the elemental energy corrections 
suggested by Lany^i for P in all phosphides 1 ^ and for Ca 
in all Ca compounds. The other elements we consider, 
Cu, Se, Zn, Cd, and CI, either have statistically insignif- 
icant corrections or, in the case of CI, are not considered 
in Ref. [nj. 

Our calculations use the projector augmented wave 
(PAW) metho d 14 ' 15 as implemented in the plane wave 
code VASP 16 . Further details of the calculation are given 
in the Supplemental Material [l7[- We use U = 6 eV 
for the Cu-d, Zn-d, and Cd-d orbitals. This value of U 
for Cu-d has been chosen in previous work (c.f. Ref. [HI] ) 
to yield agreement with the experimental band structure 
below the valence band maximum (VBM) 1 ^, thus elim- 
inating the need for post-calculation corrections to the 
VBM of Cu3PSe4^ All calculations include ionic relax- 
ation, while lattice parameters are relaxed for all pure 
compounds and elements, including the defect free host. 

The bonding character of Cu3PSe4 is evident in the 
GGA+L 7 calculation of the partial DOS, shown in Fig. [I] 
The valence bands above —7 eV and the conduction 
bands below 3.3 eV have similarities to other multinary 
copper chalcogenides. One such common property is 
that the Cu-d states are split into non-bonding e or- 
bitals and t 2 orbitals which form filled bonding and filled 
antibonding bands because of their interaction with the 
chalcogenide p orbitals^!. The antibonding band forms 
the highest valence band. Like CuInSe2, CuGaSe2, and 
Cu2ZnSnSe4^i, the conduction band has a character that 
is largely antibonding between Se-p and Mt-s, where Mt 
represents the element acting as the high valence metal 
(e.g. Sn in Cu2ZnSnSe4, P in Cu 3 PSe4). The antibond- 
ing character is inferred from the presence of a spatial 
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FIG. 2. Chemical potential domain with stable region of 
Cu3PSe4 in gray. The chosen Cu-rich growth condition is 
indicated by a circle. 



FIG. 1. (Color online) Partial density of states for the unit 
cell (Cu3PSe4)2. The large Cu-d peak rises to a maximum 
of twice the height of the plot. The vertical dashed line at 
denotes the valence band maximum. 



TABLE I. Principal axis tensor components and appropri- 
ate scalar averages for effective hole mass (units of elec- 
tron mass mo) and electronic and total dielectric constants, 
eoo and eo. The conductivity effective mass is given by 
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node between the Mt and Se atoms in the charge density 
of the lowest conduction band. Unlike materials with a 
metallic Mt, Cu3PSe4 has no valence band that is the ob- 
vious bonding counterpart. In fact, the P-s orbitals have 
nominally been filled in the P-s/Se-s bonding and anti- 
bonding bands, near —15 and —10 eV. This o~o~ bonding 
does not occur when Mt is more metallic, because of the 
larger energy difference between the atomic Mt-s level 
and the chalcogenide s level. Thus the appearance of a 
P-s/Se-p* antibond is somewhat surprising despite the 
fact that it follows the trend of other multinary copper 
chalcogenides. The bonding counterpart of the second 
conduction band, which has significant P-p/Se-p* char- 
acter, is found in the valence band near —5.7 eV. 

The calculated GGA+C/ effective hole mass and dielec- 
tric tensor components are shown in Table |TJ The dielec- 
tric tensor is calculated using density functional pertur- 
bation theory^ 2 -. The effective mass tensor, calculated 
from the band structure, has much larger components in 
the yz plane than along the x axis. Because the radius of 
a hydrogenic shallow defect state (also known as a per- 
turbed host stated) is inversely proportional to effective 
mass, this will result in shallow acceptor wavefunctions 
being greatly elongated in the x direction. The defect 
and carrier concentration calculations below do not rely 
on a value of effective mass to*, but instead use the cal- 
culated DOS directly. 



We analyze the allowed chemical potential domain for 
Cu 3 PSe 4 synthesis by calculating AH for 22 compounds 
containing Cu, P, and Se. Fig. [2] shows the results 
for several important compounds, revealing a relatively 
large stable region. To best match experimental carrier 
concentrations-, we perform the defect calculations for 
the conditions A/ip = and A/iCu = —0.11 eV (circled 
in Fig. [2]). Choosing Apcu to assume its maximum al- 
lowed value minimizes the calculated concentration of the 
shallow acceptor defect Vrj u - 

We note that it has been observed 2 ^ that under certain 
conditions Cu3PSe4 can coexist with the ionic conductor 
CuyPSeg, but this is not predicted by chemical potential 
domain analysis. This discrepancy may be due to fi- 
nite temperature effects; the low temperature a phase 24 
of Cu7PSe6 was used in calculations, while at formation 
temperature the partially disordered 7 phase would be 
present. We also note here that the error of total en- 
ergy calculations involving phosphorus can be large; a 
statistical correction of 0.6 eV per P atom is given in 
Ref. [lj| due to artifactual energy differences between 
phosphorus in reductive and neutral (elemental) environ- 
ments. This error is expected to impart uncertainty both 
to the calculated heat of formation of Cu3PSe4, which af- 
fects defect energies through its effect on Apcu, and to 
the defect supercell energies themselves, particularly for 
the high concentration Pg defect. In the latter case, 
the additional P atom is reduced by the neighboring Cu 
ions, in strong contrast to the host P atoms, which are 
oxidized by the Se neighbors. While phosphorus raises 
concern, the GGA+{7 statistical corrections 1 ^ associated 
with Cu and Se atoms are less than 0.05 eV, and our 
calculated heat of formation of Cu3Se2 is within 0.05 eV 
of experiment 10 . 

The defect analysis^ is performed initially using a 
2x2x2 (2 3 ) supercell (128 atoms). We use all finite 
size corrections^ described by Lany and Zunger— i. No 
correction is made to the valence band, while the con- 
duction band is raised so that the bandgap is increased 
from the calculated value of 0.52 eV to the experimen- 
tal value of 1.4 eV. A shallow donor correction termS 
+(-E'ff,cxp — -Eg.caic) is applied to the energies of non- 
ionized shallow donor defects. However, none of the in- 
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FIG. 3. Defect formation energies and transition energies. 
Where nonequivalent sites are calculated, the lowest energies 
for each charge state are shown. The self-consistent room 
temperature Fermi energy Ef.sc, assuming a formation tem- 
perature of 500° C, is shown by the vertical line at 0.031 eV. 
The (0/+) transition level for shallow donor Znc u has been 
raised to follow the conduction band correction. 



trinsic point defects are clearly shallow donors, and thus 
this correction is applied only for the extrinsic donors 
considered: Ca, Cd, and Zn on a Cu site, and CI on a Se 
site. As a result, none of the intrinsic defect transition 
energies, which are given with respect to the VBM, have 
been adjusted with the widening of the gap. 

Formation energies and transition levels for the lower 
energy intrinsic defects and the lowest energy extrinsic 
defect are shown in Fig. [3] The acceptors Vc u and Ps e 
both pin the Fermi energy below mid-gap, preventing 
Cu3PSe4 from being n-doped near thermal equilibrium. 
The formation energy of the neutral defect V^ u is calcu- 
lated to be 0.50 eV, with a (-/0) transition level of 0.05 to 
0.06 eV, depending on the Cu site. The formation energy 
of Pg varies with site from 0.47 eV to 0.50 eV, with the 
(-/0) transition levels varying from 0.08 eV to 0.17 eV. 

We use the formation temperature of 500°C (approxi- 
mately the temperature used in recent pellet and single 
crystal experiments^) to calculate^ the self-consistent 
formation temperature defect population and Fermi level. 
The resulting total defect concentrations (irrespective of 
charge state) are 4.1 x 10 19 cm~ 3 for Vc u and 4.6 x 10 19 
cm -3 for Pg . We hx the total defect concentrations and 
recalculate the self-consistent populations of the various 
charge states at room temperature (300 K). The room 
temperature self-consistent Fermi level is 0.031 eV above 
the VBM with a hole concentration of p = 8 x 10 18 cm~ 3 . 
The contribution of V^ u to the hole density is over five 
times that of Pg c , indicating that Vq u is electronically 
the most important defect type. 

If zinc is present during formation, the maximum ZnJ u 
concentration is approximately 5 x 10 18 cm~ 3 , and the 
net room temperature hole density is lowered slightly to 



TABLE II. Calculated site-averaged formation energies for 
Vcu defects and predicted versus experimental hole concen- 
trations. The error of the GGA+U method is seen to be 
smaller than the more standard^ method of using GGA in- 
cluding VBM corrections. 
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p = 6 x 10 18 cm~ 3 . The other potential donor dopants 
considered have greater formation energies and can be 
neglected. 

The Vcu defect is very weakly localized and we have 
recalculated its two charge states on one Cu site using 
a 4 3 (1024 atom) supercell. Even for this supercell size, 
the defect wavefunction is not localized within the super- 
cell in the x direction (the low effective mass direction). 
We obtain AH = 0.53 eV for the neutral defect, with 
a (-/0) charge transition level of 0.04 eV. We note that 
the hydrogen-like approximation using the conductivity 
effective mass m* ond = 0.27 mo yields a comparable bind- 
ing energy of 0.02 eV. Assigning the large supercell data 
to all Cu sites, in combination with the previous Ps c data, 
yields an insignificantly modified hole density p = 9 x 10 18 
cm~ 3 . 

The Psc defect state is substantially localized within 
the smaller 2 3 supercell and has P-p character on the 
defect site and Cu-d character on the nearest neighbors. 
The degree of localization allows us to apply the defect 
image charge correction of Refs. [1] using the neutral de- 
fect potential as the reference potential^. The resulting 
correction (0.08 eV) agrees well with the corresponding 
correction (image charge^ + alignment = 0.07 eV) ac- 
cording to Refs. [1 and[3|. 

The defect analysis performed here agrees qualitatively 
with recent experimental results. Our calculated Cu:P 
ratio of 2.97 is consistent with the value 2.92 ± 0.06 mea- 
sured for single crystals 8 . We predict a large hole con- 
centration of p = 8 x 10 18 cm~ 3 , about one order of 
magnitude larger than the value 6 x 10 17 cm -3 obtained 
by Hall and Seebeck measurements on pressed, sintered 
pellets 8 -. 

We compare the GGA+U defect calculations described 
above with standard GGA defect calculations followed 
by application of a GGA+J7 correction&i 8 - (—0.34 eV) 
to the VBM. The GGA defect calculations include all 
types of corrections applied to the GGA+U calculations 
and include a GGA determination of the maximum al- 
lowed copper chemical potential A/icu (—0.06 eV). As 
shown in Table HH the more standard "GGA + VBM 
correction" procedure changes the formation energies of 
Vq u and Vq u (evaluated at maximum A/icu and min- 
imum Ep ) by about —0.16 eV, causing a significantly 
larger overestimation of p. This comparison indicates 
that the superiority of GGA+U to GGA in bulk total 
energy calculation a 10 ' 11 is also extended to defect calcu- 
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lations. 

It is instructive to consider further the implications 
of the available experimental results 8 -. We estimate the 
amount by which AH values would have to change in 
order to bring the calculated hole concentration p closer 
to the value measured for polycrystalline pellets 8 -. If one 
assumes that the calculated transition energy of Vc u is 
not underestimated, the experiments of Ref. [8j indicate 
that the formation energy of Vc u must increase, while 
the transition level of Ps c increases and the formation 
energy of Ps e decreases. The adjustment to the Vo en- 
ergy must be significant to recover the measured p. For 
example, increasing the formation energy of Vc u defects 
by 0.35 eV while applying changes of —0.05 and 0.05 eV 
to the neutral and charged Pse defects respectively yields 
p = 7 x 10 17 cm~ 3 and a Cu:P ratio of 2.96. Such large 
changes to the Vq u formation energies are not easily ex- 
plained by calculational errors associated primarily with 
phosphorus. 

An alternative possibility is that the GGA+U calcu- 
lated VBM is too high by a modest amount, and that 
the apparent shallow character of the Vc u defect is an 
artifact of this band misplacement. For example, ap- 
plying a valence band correction AEy — —0.1 eV and 
choosing not to apply the shallow acceptor corrections to 
the neutral defects (that is, using a strictly "band edge 
only" approach) yields the much lower hole concentra- 
tion p=l.lx 10 18 cm~ 3 with an only slightly increased 
Cu:P ratio (2.976). 

The experimental data suggests an increase in neutral 
P concentration, and possibly the presence of a low en- 
ergy donor defect involving extra P atoms, which could 
lower p by compensating the Vc u acceptors and thus 
avoid the need to raise AiJ(Vcu)- We therefore have 
examined^, at lower accuracy and without finite size cor- 
rections, several configurations of neutral P intcrstitials 
Ip as well as a number of neutral two-defect complexes 



of types Pse-Pse, PseTp, and Ip-Ip. The lowest complex 
formation energy that is less than the sum of its two parts 
is 1.5 eV for a Pse-Ip complex. The Ip formation energies 
range from 1.8 to 2.4 eV. Additionally, we have examined 
several of these defects in the +1 (active donor) charge 
state but have not obtained energies low enough to com- 
pete with the calculated AH(V Cu ) = 0.56 eV. Thus no P 
interstitial or P-rich defect complex in the energy range 
of interest has been found. 

In conclusion, we have performed a set of GGA+C/ 
defect calculations on Cu3PSe4, a p-type semiconduc- 
tor with a direct bandgap of 1.4 eV. We compare our 
methods against standard GGA, larger supercells, and 
alternative correction methods. We predict that the Vc u 
defect is mostly responsible for the large, experimentally 
observed intrinsic hole concentration p, with some con- 
tribution from Pse- Both of these defects pin the Fermi 
level below mid-gap, so that n-doping is prohibited near 
thermal equilibrium. Both defects also contribute to the 
observed non-stoichiometric Cu:P ratio. Our calculation 
overestimates p by about one order of magnitude. How- 
ever, the present GGA+f/ method is shown to be more 
accurate than more standard GGA calculations with va- 
lence band corrections. Doping with Zn is calculated to 
have a small but noticeable effect on p. Because of the 
apparent uncertainty in the calculations however, this 
analysis does not rule out the possibility that Zn doping 
could significantly reduce p. 
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The two-thirds monopole term-' is used for the electrostatic 
image corrections. The potential alignment AVpa for both 
neutral and charged defects is calculated using an average 
of core potentials excluding the defect site, and in some 
cases nearest neighbor sites, as described in Ref. 0]. This 
alignment is applied as a correction (4-gAVpA) for charged 
defect energies, and is also used to adjust the values of the 
host VBM and conduction band maximum (CBM) that 
are used in the band filling correction described in Ref. [(| • 



